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Abstract 

We compare the performance of extremal optimization (EO), flat- 
histogram and equal-hit algorithms for finding spin-glass ground states. 
The first-passage-times to a ground state are computed. At optimal pa- 
rameter of r = 1.15, EO outperforms other methods for small system sizes, 
but equal-hit algorithm is competitive to EO, particularly for large sys- 
tems. Flat-histogram and equal-hit algorithms offer additional advantage 
that they can be used for equilibrium thermodynamic calculations. We 
also propose a method to turn EO into a useful algorithm for equilibrium 
calculations. 

Keywords: extremal optimization, flat-histogram algorithm, equal-hit 
algorithm, spin-glass model, ground state. 

1 Introduction 

Optimization with methods motivated from real physical processes is an active 
field of research. Simulated annealing and genetic algorithm || are two 
well-known examples. In particular, there have been a large variety of methods 
proposed to find spin-glass ground states || ||, [|, 0, JtL |L [|, ffo|. Recently, 
Boettcher and Percus [jll| [lj] introduced 'extremal optimization' (EO) inspired 
by models of self-organized criticality p3| , which gave impressive performance. 

Most of the heuristic optimization methods (including simulated annealing 
and genetic algorithm) are designed to find ground states only, thus it is not 
possible to give correct thermodynamics from a simulation. On the other hand, 
multi-canonical ensemble simulation l/fc-sampling |l5| ], parallel tempering 
[ fl6| , and recent flat-histogram dynamics [[l7| are constructed for equilibrium 
thermodynamics, but can also be used as methods for optimization. A study 



of optimization by flat-histogram algorithm on the two-dimensional spin-glass 
model is carried out in rcf . ]f8| . Unlike simulated annealing and other heuristic 
methods, we note that these methods do not have any adjustable parameters. 
It is useful to know the efficiencies of this second class of methods when used as 
an optimization tool. 

In this paper, we make a comparative study of the extremal optimization and 
flat-histogram/equal- hit dynamics. We compare four algorithms: EO at r = 1 
with a continuous approximation in the probability of choosing a site, original 
EO at optimal value of r = 1.15, single-spin-flip flat-histogram dynamics, and 
equal- hit algorithm with iV-fold way. It is found that EO at the value r = 1.15 
is very good for both two- and three-dimensional Ising spin glasses. The equal- 
hit algorithm with TV-fold way is also competitive. For large systems, equal-hit 
appears even slightly better than EO. It is useful to have the efficiency of EO 
but still give equilibrium results. To this end, we introduce a rejection step in 
EO, thus turning EO into an equilibrium simulation method. 



2 Single-spin-flip algorithms 

In the following, we specialize our discussion in the context of spin models, and 
particularly the spin-glass model |19| . The energy function is defined by 



E ( a ) = ~Y1 J y cr ' cr J' (!) 

{id) 



where the spin cr^ takes on value +1 or —1 with i varying over a hypercubic 
lattice in d dimensions. The coupling constant J,j for each nearest neighbor 
pair takes on a random value of +J and — J with equal probability. We 
impose a constraint Y^,/ij) Jij = 0- The spin glass is known to be one of the 
hardest problems [^0| to find the state a that minimizes E{a). 

A single-spin-flip with rejection is described by a Markov chain transition 
matrix of the form 

W(a -> a') = S N (a, <j')^a(a - a'), a ? a', (2) 

where <5jv(c, a') = 1 if a' is obtained from a by a single spin flip, and otherwise. 
The factor 1/N represents the random selection of a spin, where N (= L d ) is the 
number of spins. a(a — > a') is the flip rate. If we choose a(cr — > a') according 
to Metropolis rate, 



. / f(E(a')) 
a((j — > a ) = mm 1, — 7 — , V 



(3) 



we can realize equilibrium distribution with the probability of states distributed 
according to /(£J(cr)J. Some choices are: 

exp(— E / (ksT)) , canonical ensemble; 
f(E) = I l/n>{E), multicanonical ensemble; (4) 

1/ Jf n(E')dE', l/k sampling, 
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where T is temperature, ks is Boltzmann constant, and n(E) is density of states 
at energy E. 

Arbitrary choice of the flip rate a(a — ► a') in general would not give one 
important property of the equilibrium systems, i.e., the microcanonical property 
that the probability distribution of the configuration <r is a function of energy 
E only. For example, the original broad histogram rate plj 

a(a^*)=rmn(l,^- ( - r ), (5) 

and the random walk rate of Berg p2[ do not have microcanonical property, 
where Nk{<r) is the number of possible moves of class k in the state a; we 
associate a class for each site i with a number from to Z (= 2c?) by a scaled 
energy change k = ((E(a') - E{a))/J + 2Z)/4, i.e., 
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yCTiO-j + 1). (6) 



7Vfe(cr) is the number of such sites having a class k. 

The single-spin-flip version with rejection can be turned into a rejection- free 
iV-fold way p3f simulation where a class is chosen with probability 

Pk = ^T> A(a)=J2a k N k (a), (7) 

where a k is a(cr — > cr') associated with an energy change indicated by class k. 
A spin in that class is picked up at random, and the flip is always accepted. 
Thermodynamic quantities need to be weighted by a factor 1/A(a). 

EO Q |l2| is somewhat related to iV-fold way in the sense that a class 
is chosen with some probability Pk, and a spin in that class is picked up and 
always flipped. The EO algorithm can be stated as follows: we classify the site 
by its 'fitness' k. There are Z + 1 possible values for k. In the general EO, the 
sites are sorted according to the fitness k. Since there are only a small number 
of possible values in the ±J spin- glass model, the sorting is not necessary. We 
simply make a list of sites in each category. We pick a class according to the 
probability P k , and then choose a spin in that class and flip with probability 
one. The corresponding transition matrix is then 

W(a^a') =S N (a,a')P k ^, a + & '. (8) 

The original choice of EO is to take 

Na+N 1 + ---+N k - 1 <i<N + N 1 ^ hW fc 

with t being a parameter of the algorithm. We define the standard EO to be 
a continuous approximation to the above sum at r = 1 with the analytical 
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expression by 



1 ,„i±St 



TV, 



^MNTT^T^W,- (10) 

The number 1 in the numerator and denominator are introduced somewhat 
arbitrarily to avoid divergence when the sum of Nk is zero. An optimized EO 
will be the discrete version, Eq. (0), with r that gives best performance; we use 
r = 1.15 as recommended in ref [|12| . To realize the power-law distribution, we 
generate an integer i = [^ 1 /( 1_T )J , where £ is a uniformly distributed random 
number between and 1, and pick a corresponding site i ordered by the class. 
We have used [• • -J for the floor function. 



3 Comparison of EO with flat-histogram and 
equal-hit algorithms 

The flat-histogram algorithm jl7], ^| is a special choice of the flip rate 



FH7 r\ ■ I i ( N Z~k{<j'))E' i , .... , 

a =™(^ {Nk{(j))E ), (ID 

where the angular brackets with subscript E denote a microcanonical average of 
the quantity TVfc(cr) at energy E; the starting state a has energy E and the final 
state a' has energy E' . This particular choice of the rate gives a flat distribution 
for the energy histogram, H(E) = n(E)f(E) — const, n(E) is density of states. 
This is one way to realize multicanonical ensemble. 

In the TV-fold way equal-hit algorithm |2^] , we perform the usual TV-fold- 
way move (thus rejection-free) which is constructed from the following single- 
spin-flip rate: 

,EQ/_ . _^_„ s _ A (A) E {N Z -k(v'))i 



a^(a a') = min 1, v ^ 7 ~ / " ■ (12) 
V (A) E >(N k (a)) E J 

We note that (A}^ 1 = {1/A)n, where (• • -)n is average over the TV-fold way 
samples. In equal-hit algorithm, it is guaranteed by construction that we change 
state in every move, and the distribution of the visits to different energies is flat. 

Since the microcanonical averages used in the flip rates are not known before 
the simulation, we use running average to replace the exact microcanonical 
average. It appears that this is a valid approximation and should converge to 
the correct values for sufficiently long runs. For a truly exact algorithm (in 
the sense of realizing microcanonical property), it is sufficient with a two-pass 
simulation. The first pass uses a running average; in the second pass, we use a 
multicanonical rate determined from the first pass. 

Many different criteria are used to measure the effectiveness of an optimiza- 
tion algorithm, such as the fraction of cases for which ground states are found 
in a set of runs. The first-passage-time, the time in units of Monte Carlo sweeps 
that a ground state is found, starting with similar random configurations, should 
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Figure 1 : The average first-passage time t g in units of sweeps to find a ground 
state for four algorithms: standard EO (triangles), optimal EO (circles), flat- 
histogram (diamonds), and equal-hit (squares) for the two-dimensional spin- 
glass model. Over 10 3 realization of random coupling samples are used for 
averaging for each algorithm and size. 



be a good measure of the algorithms' efficiency. We consider sample average of 
the first-passage-time, although the distribution of it is also very useful. The 
computer CPU time is another useful criterion when comparing algorithms of 
very different types. 

In the flat-histogram and equal- hit algorithms, we can sample positive as 
well as negative energies uniformly. In this study, we have restricted to the 
negative energy part, where moves to E > region are rejected. We compute 
the average time (first-passage time) for each lattice size and given algorithm 
in units of sweeps (N = L d basic moves) to find a ground state, starting from a 
random configuration of equal probability of spin up and spin down. For two- 
dimensional ± J spin glass, we determine the first-passage time t g by comparing 
the current value of energy with the exact value of ground state energy, obtained 
from the Spin Glass Server p6| . Thus the results are unbiased. The average 
first-passage time t g for the two-dimensional Ising spin-glass model is shown in 
Fig. |l|. Over 10 3 realization of random coupling samples are used for averaging 
for each algorithm and size. 

Since the ground state energies are usually not known in three dimensions, 
we consider instead the time for finding the lowest energy for a fixed amount of 



5 



10° 



10° 



3D spin glass 



EO,tau=1 



c» 10 



10 2 



10° 



/ 



/ / 



Flat-Histog^rfi 



• Equal-Hit 



EO, tau=1.15 



6 8 10 12 

L 



Figure 2: The average limiting time t g to find a ground state for the three- 
dimensional Ising spin-glass model. The meanings of the symbols are the same 
as those of Figure 1 . 

sweeps t, averaged over the coupling constants with the constraint ^ Jy = 0. 
For any fixed running length t, results obtained are only a lower bound for 
t g . We consider run lengths of 10 4 , fO 5 , 10 6 , etc, until the first-passage time 
converges for large t. This limiting time t g is reported for the three-dimensional 
Ising spin-glass model in Fig. |2[ 

To compare the efficiency of the four algorithms, the actual CPU times are 
also an important factor. For our implementation, it turns out that the opti- 
mized EO, standard EO, and iV-fold way equal-hit all have about the same speed 
at 6 microsecond per spin flip on a 700 MHz Pentium, while the single-spin-flip 
flat-histogram algorithm takes 3 microsecond. There are several important fea- 
tures in this comparison, see Fig. [|. All of them have a first-passage time that 
grows rapidly with sizes. With the exception of the standard EO, they nearly 
have the same slope of about 6 on a double logarithmic scale. It is also interest- 
ing to compare the first-passage time with that of equilibrium tunneling time 
reported in ref. |l8|, |24| . EO gives excellent performance for small to moderate 
size systems. However, for large sizes, equal-hit is as good as EO, or even bet- 
ter. Flat-histogram is worse by some constant factor. On the other hand, the 
performance of the standard EO at r = 1 is rather poor. This shows that the 
results of EO is rather sensitive to the value of t. 

Another very interesting aspect of Fig. |l| is that the curves all look linear in 
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L = Q, MCS = 10 e , sample = 1024 


EO 

optimal EO 

Flat-Histogram 

Equal-Hit 


(6.95 ±0.65) x 10 4 
(1.36 ±0.15) x 10 3 
(1.46 ±0.13) x 10 4 
(1.92 ±0.20) x 10 3 


E g 

-1.7713 ±0.0012 
-1.7715 ±0.0006 
-1.7715 ±0.0006 
-1.7716 ±0.0006 


L = 10, MCS = 10 b , sample = 256 


EO 

optimal EO 

Flat-Histogram 

Equal-Hit 


h 

(2.98 ±0.17) x 10 5 
(1.04 ±0.11) x 10 5 
(1.93 ±0.15) x 10 5 
(1.27 ±0.17) x 10 5 


E 9 

-1.7721 ±0.0008 
-1.7809 ±0.0008 
-1.7802 ±0.0007 
-1.7815 ±0.0010 



Table 1: The average time t g and lowest energy E g obtained for three- 
dimensional Ising spin-glass. 

the semi-logarithmic scale. This implies that t g ~ exp(cL) for some constant 
c ~ 1, not a power law in L. Thus, all of the algorithms are asymptotically 
inefficient. It would be very interesting if this numerical observation can be 
supported by some argument. Similar results for the three-dimensional Ising 
spin glass is presented in Fig. ||. 

In Table [l], we report some typical data for the average first-passage time 
t g , energy per site, length of the run, and number of random samples for four 
algorithms for the three-dimensional spin glasses. Since we use the same set of 
samples with the four algorithms, a lower energy indicates a better performance. 
The data show that equal-hit is comparable to EO at optimal r. 

4 Turning EO into an equilibrium algorithm 

The flat-histogram and equal-hit algorithms can be used for equilibrium simu- 
lation. With the help of counting the number of potential moves, N k , a basic 
requirement for obtaining equilibrium property of the simulated model is the 
microcanonical property. Using the broad histogram equation j2?j, [2^] , 

n(E){N k (a)) E =n{E'){N z _ k (a')) E ,, k = {{E 1 - E)/J + 2Z)/4, (13) 

we can obtain the density of states n(E) of energy E, thus the equilibrium 
thermodynamic quantities, including free energy. 

Unfortunately, the microcanonical property that the probability distribution 
P(cr) of the configurations is a function of energy E only is strongly violated 
in EO. The probabilities of the ground states cluster into groups, rather than 
uniformly distributed. Numerical tests show that P(cr) is a function of E, N k) as 
well as additional unknown parameters. To correct this problem, we introduce 
a rejection step in the EO algorithm, as follows: 

W(a^a / )=S N (a,a / )P k ^-a(a^a'). (14) 



7 



The acceptance rate a is determined by imposing a detailed balance with an 
unknown probability distribution f(E(a)), 

f(E(a))W(a - a') = f(E(a'))W(a' -» a). (15) 

This gives an equation for the rate a: 

f{E)P k -±—a{o - O = f(E')P' z _ k 1 a(a' -> a). (16) 

The prime on P' indicates that it is a set of P values calculated from the state 
a' . A solution to this equation is a Metropolis-type choice: 

^ ° ] = mm I'' f(E)P k /N k (a) ) • < 17 > 

To implement this, we need a two-pass simulation. The first pass determines 
the function f{E). The procedure is by no means unique. Here, we collect 
histogram of energy H(E) as well as statistics for (N k )E from an incorrect 
simulate of the original EO. Then we determine an approximate density of 
states with the help of the broad histogram equation, Eq. (|l^) . The function / 
is computed from f(E) = H(E)/n(E). The EO with rejection is implemented in 
the second pass. The above procedure should be applicable for any model and 
any optimization algorithms that has a 'steady state'. If the microcanonical 
property is only slightly violated, it will give a correct equilibrium algorithm 
with a nearly equal to 1. Thus, we hope to have a method that is efficient for 
optimization, and yet at the same time, give correct equilibrium statistics. 

Indeed, with the above method the microcanonical property is restored. Due 
to the rejection step, the dynamics is slightly changed. A consequence is that 
the histogram in the second pass shifted towards high energy side, thus the 
efficiency of the original EO is lost. 



5 Conclusion 

From this study, we show that equal-hit algorithm is an excellent candidate for 
ground state search. At the same time, it also offers the possibility for equi- 
librium calculations, such as the computation of the ground state entropy. We 
also show how optimization algorithms like EO can be turned into equilibrium 
algorithms by introducing a rejection step. All the algorithms studied here give 
rather rapid increase of t g with sizes, thus it is important and challenging to 
find algorithms that reduce this growth. Perhaps, algorithms based on single- 
spin-flip have their fundamental limitations. 
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